Deforestation Risk

Abstract

A workflow for deriving jurisidictional risk-allocated deforestation mapping compliant with Verra’s VMD0055 (V1.1) module and the VM0048 (V1.0) consolidated methodology. This includes baseline estimates of forest change and jurisdictional allocated deforestation.

Keywords

REDD+, VCS, Verra, Carbon verification, Jurisdictional

1. Introduction

This assessment develops a spatially explicit deforestation risk map and allocation framework for Bong County, Liberia, following Verra’s VM0048 methodology requirements. The analysis integrates the baseline emissions estimates developed in the previous workflow with spatial risk modeling to allocate jurisdictional deforestation across the landscape based on empirically-derived risk factors.

This analysis maintains consistency with the baseline emissions assessment by focusing exclusively on Bong County (8,772 km²), located in west-central Liberia. The county encompasses diverse forest-agricultural landscapes within the Upper Guinea Forest biodiversity hotspot, representing typical West African deforestation pressures including:

• Agricultural expansion (primarily rice and cassava)
• Small-scale logging and charcoal production
• Road development and settlement expansion
• Mining activities (iron ore and artisan mining)

2. Method

Four-step process:

1. Covariate Development: Process infrastructure, demographic, and biophysical risk data layers
2. Statistical Evaluation: Check covariate magnitude to historical deforestation
3. Model Development: Empirically weighted risk model generation
4. Deforestation Allocation: Spatial distribution of baseline deforestation rates based on risk surfaces

# Import & reproject to OSM grid systems
aoi_country = geodata::gadm(country="LBR", level=0, path="./assets/AOI") |> sf::st_as_sf() |> sf::st_transform("EPSG:3857")
aoi_states  = geodata::gadm(country="LBR", level=1, path="./assets/AOI") |> sf::st_as_sf() |>
  dplyr::filter(NAME_1 == "Bong") |> sf::st_transform("EPSG:3857")
crs_master = sf::st_crs(aoi_country)
  
# Total extent
area_km2 = as.numeric(sf::st_area(aoi_states)) / 1e6
cat("Study Area: Bong County, Liberia\n")

NA Study Area: Bong County, Liberia

cat("Study Area:", round(area_km2, 0), "km²\n")

NA Study Area: 8565 km²

# Visualize AOI
tmap::tmap_mode("view")
tmap::tm_shape(aoi_country) + tmap::tm_borders(col="purple", lwd=2) +
  tmap::tm_shape(aoi_states) + tmap::tm_borders(col="red", lwd=1) +
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
  tmap::tm_title("AOI Jurisdictional Boundaries", size=.8) + 
  tmap::tm_basemap("OpenTopoMap") 

Area check

In Liberia, the official definition of forest land is provided by the Forestry Development Authority (Government of Liberia 2019), including areas of land that meet the following criteria:

• Canopy cover of minimum 30%;
• Canopy height of minimum 5m or the capacity to reach it;
• Covering a minimum of 1 hectare of land.

forest_2024=terra::rast("./assets/LULC/outputs/forest_2024.tif")
forest_2024_poly=terra::as.polygons(forest_2024) |> sf::st_as_sf()
forest_2024_poly$area_ha = round(as.numeric(sf::st_area(forest_2024_poly) * 0.0001, 4))
forest_2024_poly |> sf::st_drop_geometry() |> janitor::adorn_totals() 

label	area_ha
0	33214
1	805113
Total	838327

slivers = forest_2024_poly |> dplyr::filter(as.numeric(area_ha) < 1)  
slivers 

label	geometry	area_ha

# no artefacts found??

---

Spatial Covariates

The following spatial covariates were processed as potential drivers of deforestation risk. Covariates were merged between sociodemographic and geographic datasets surrounding the project area and national level datasets beyond the project area in order to enable jurisdictions analysis.

Distance Covariates

# Distance to forest edge
forest_2024 = terra::rast("./assets/LULC/outputs/forest_2024.tif")
forest <- terra::classify(forest_2024, cbind(c(1,0), c(1,NA)))
#nonForest <- terra::classify(forest_2024, cbind(c(1,0), c(NA,1)))
#writeRaster(forest, "./assets/LULC/outputs/forest.tif", overwrite=T)
#writeRaster(nonForest, "./assets/LULC/outputs/nonForest.tif", overwrite=T)
forest_edge = terra::boundaries(forest, inner=FALSE)
forest_edge_distance = terra::distance(forest_edge)
forest_edge_distance = forest_edge_distance / 1000
names(forest_edge_distance) = "edge_distance_km"
forest_edge_distance

# Distance to road # Check OSM data & derive bbox window
osmdata::available_features()  
osmdata::available_tags("route")
# "highway", "public_transport", 
# "railway", "route", 
# "side_road", "tracks"
osm_bbox = osmdata::getbb("Liberia")
roads_osm = osm_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "route") |>
  osmdata::osmdata_sf()

# guardrail
if(!is.null(roads_osm$osm_lines)) {
  roads_bong = roads_osm$osm_lines |> 
    sf::st_intersection(aoi_states) |>
    sf::st_transform("EPSG:3857")

# Hospitals
buildings_hospital = aoi_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "amenity", value = "hospital") |>
  osmdata::osmdata_sf() # |> st_transform(crs_master)
roads_rast <- rasterize(vect(transport), template, field=1, background=0)
roads_rast <- subst(roads_rast, from=0, to=NA)
distance_to_roads <- distance(roads_rast)
writeRaster(distance_to_roads, "./data/LULC/distance_to_roads.tif", overwrite=T)

Error in parse(text = input): <text>:40:0: unexpected end of input
38: 
39: 
   ^

Built Environment

# Check OSM data & derive bbox window
osmdata::available_features()  

NA   [1] "4wd_only"                    "abandoned"                  
NA   [3] "abutters"                    "access"                     
NA   [5] "addr"                        "addr:*"                     
NA   [7] "addr:city"                   "addr:conscriptionnumber"    
NA   [9] "addr:country"                "addr:county"                
NA  [11] "addr:district"               "addr:flats"                 
NA  [13] "addr:full"                   "addr:hamlet"                
NA  [15] "addr:housename"              "addr:housenumber"           
NA  [17] "addr:inclusion"              "addr:interpolation"         
NA  [19] "addr:place"                  "addr:postbox"               
NA  [21] "addr:postcode"               "addr:province"              
NA  [23] "addr:state"                  "addr:street"                
NA  [25] "addr:subdistrict"            "addr:suburb"                
NA  [27] "addr:unit"                   "admin_level"                
NA  [29] "aeroway"                     "agricultural"               
NA  [31] "alcohol"                     "alt_name"                   
NA  [33] "amenity"                     "area"                       
NA  [35] "atv"                         "backward"                   
NA  [37] "barrier"                     "basin"                      
NA  [39] "bdouble"                     "bicycle"                    
NA  [41] "bicycle_road"                "biergarten"                 
NA  [43] "boat"                        "border_type"                
NA  [45] "boundary"                    "brand"                      
NA  [47] "bridge"                      "building"                   
NA  [49] "building:colour"             "building:fireproof"         
NA  [51] "building:flats"              "building:levels"            
NA  [53] "building:material"           "building:min_level"         
NA  [55] "building:part"               "building:soft_storey"       
NA  [57] "bus"                         "bus_bay"                    
NA  [59] "bus:lanes"                   "busway"                     
NA  [61] "capacity"                    "carriage"                   
NA  [63] "castle_type"                 "change"                     
NA  [65] "charge"                      "clothes"                    
NA  [67] "construction"                "construction_date"          
NA  [69] "construction#Railways"       "covered"                    
NA  [71] "craft"                       "crossing"                   
NA  [73] "crossing:island"             "cuisine"                    
NA  [75] "cutting"                     "cycle_rickshaw"             
NA  [77] "cycleway"                    "cycleway:left"              
NA  [79] "cycleway:left:oneway"        "cycleway:right"             
NA  [81] "cycleway:right:oneway"       "denomination"               
NA  [83] "destination"                 "diet:*"                     
NA  [85] "direction"                   "dispensing"                 
NA  [87] "disused"                     "dog"                        
NA  [89] "drinking_water"              "drinking_water:legal"       
NA  [91] "drive_in"                    "drive_through"              
NA  [93] "ele"                         "electric_bicycle"           
NA  [95] "electrified"                 "embankment"                 
NA  [97] "embedded_rails"              "emergency"                  
NA  [99] "end_date"                    "energy_class"               
NA [101] "entrance"                    "est_width"                  
NA [103] "fee"                         "female"                     
NA [105] "fire_object:type"            "fire_operator"              
NA [107] "fire_rank"                   "food"                       
NA [109] "foot"                        "footway"                    
NA [111] "ford"                        "forestry"                   
NA [113] "forward"                     "frequency"                  
NA [115] "frontage_road"               "fuel"                       
NA [117] "full_name"                   "gauge"                      
NA [119] "gender_segregated"           "golf_cart"                  
NA [121] "goods"                       "gutter"                     
NA [123] "hand_cart"                   "hazard"                     
NA [125] "hazmat"                      "healthcare"                 
NA [127] "healthcare:counselling"      "healthcare:speciality"      
NA [129] "height"                      "hgv"                        
NA [131] "highway"                     "historic"                   
NA [133] "horse"                       "hot_water"                  
NA [135] "hov"                         "ice_road"                   
NA [137] "incline"                     "industrial"                 
NA [139] "inline_skates"               "inscription"                
NA [141] "int_name"                    "internet_access"            
NA [143] "junction"                    "kerb"                       
NA [145] "landuse"                     "lane_markings"              
NA [147] "lanes"                       "lanes:bus"                  
NA [149] "lanes:psv"                   "layer"                      
NA [151] "leaf_cycle"                  "leaf_type"                  
NA [153] "leisure"                     "lhv"                        
NA [155] "lit"                         "loc_name"                   
NA [157] "location"                    "male"                       
NA [159] "man_made"                    "max_age"                    
NA [161] "max_level"                   "maxaxleload"                
NA [163] "maxheight"                   "maxlength"                  
NA [165] "maxspeed"                    "maxstay"                    
NA [167] "maxweight"                   "maxwidth"                   
NA [169] "military"                    "min_age"                    
NA [171] "min_level"                   "minspeed"                   
NA [173] "mofa"                        "moped"                      
NA [175] "motor_vehicle"               "motorboat"                  
NA [177] "motorcar"                    "motorcycle"                 
NA [179] "motorroad"                   "mountain_pass"              
NA [181] "mtb:description"             "mtb:scale"                  
NA [183] "name"                        "name_1"                     
NA [185] "name_2"                      "name:left"                  
NA [187] "name:right"                  "narrow"                     
NA [189] "nat_name"                    "natural"                    
NA [191] "nickname"                    "noexit"                     
NA [193] "non_existent_levels"         "nudism"                     
NA [195] "office"                      "official_name"              
NA [197] "old_name"                    "oneway"                     
NA [199] "oneway:bicycle"              "oneway:bus"                 
NA [201] "openfire"                    "opening_hours"              
NA [203] "opening_hours:drive_through" "operator"                   
NA [205] "orientation"                 "oven"                       
NA [207] "overtaking"                  "parking"                    
NA [209] "parking:condition"           "parking:lane"               
NA [211] "passenger_lines"             "passing_places"             
NA [213] "place"                       "power"                      
NA [215] "power_supply"                "priority"                   
NA [217] "priority_road"               "produce"                    
NA [219] "proposed"                    "proposed:name"              
NA [221] "protected_area"              "psv"                        
NA [223] "psv:lanes"                   "public_transport"           
NA [225] "railway"                     "railway:preserved"          
NA [227] "railway:track_ref"           "recycling_type"             
NA [229] "ref"                         "ref_name"                   
NA [231] "reg_name"                    "religion"                   
NA [233] "religious_level"             "rental"                     
NA [235] "residential"                 "roadtrain"                  
NA [237] "route"                       "sac_scale"                  
NA [239] "sauna"                       "service"                    
NA [241] "service_times"               "shelter_type"               
NA [243] "shop"                        "short_name"                 
NA [245] "shoulder"                    "shower"                     
NA [247] "side_road"                   "sidewalk"                   
NA [249] "site"                        "ski"                        
NA [251] "smoking"                     "smoothness"                 
NA [253] "social_facility"             "sorting_name"               
NA [255] "speed_pedelec"               "sport"                      
NA [257] "start_date"                  "step_count"                 
NA [259] "substation"                  "surface"                    
NA [261] "tactile_paving"              "tank"                       
NA [263] "taxi"                        "tidal"                      
NA [265] "toilets"                     "toilets:wheelchair"         
NA [267] "toll"                        "topless"                    
NA [269] "tourism"                     "tourist_bus"                
NA [271] "tracks"                      "tracktype"                  
NA [273] "traffic_calming"             "traffic_sign"               
NA [275] "trail_visibility"            "trailblazed"                
NA [277] "trailblazed:visibility"      "trailer"                    
NA [279] "tunnel"                      "turn"                       
NA [281] "type"                        "unisex"                     
NA [283] "usage"                       "vehicle"                    
NA [285] "vending"                     "voltage"                    
NA [287] "water"                       "wheelchair"                 
NA [289] "wholesale"                   "width"                      
NA [291] "winter_road"                 "wood"

osmdata::available_tags("place")

Key	Value
place	allotments
place	archipelago
place	borough
place	city
place	city_block
place	continent
place	country
place	county
place	district
place	farm
place	hamlet
place	island
place	islet
place	isolated_dwelling
place	locality
place	municipality
place	neighbourhood
place	ocean
place	plot
place	polder
place	province
place	quarter
place	region
place	sea
place	square
place	state
place	subdistrict
place	suburb
place	town
place	village

aoi_bbox = osmdata::getbb("Liberia")

# Hospitals
buildings_hospital = aoi_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "amenity", value = "hospital") |>
  osmdata::osmdata_sf() # |> st_transform(crs_master)

# Housing
buildings_residential = aoi_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "building", value = "residential") |>
  osmdata::osmdata_sf() # |> st_transform(crs_master)

# Religion
buildings_workship = aoi_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "building", value = "religious") |>
  osmdata::osmdata_sf() 

# Towns 
townnames = aoi_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "place", value = "town") |>
  osmdata::osmdata_sf()

# Villages
villages = aoi_bbox |> osmdata::opq() |>
  osmdata::add_osm_feature(key = "place", value = "village") |>
  osmdata::osmdata_sf()

# Merging built collections
places_merged = townnames$osm_points |> bind_rows(villages$osm_points) |>
  group_by(across(-geometry)) |> sf::st_cast("POINT") |>
  summarise(geometry = st_union(geometry), .groups = "drop")
sf::st_write(places_merged, "./03_Spatial_Data/POP/places_merged.shp", delete_dsn=T) 

NA Deleting source `./03_Spatial_Data/POP/places_merged.shp' failed
NA Writing layer `places_merged' to data source 
NA   `./03_Spatial_Data/POP/places_merged.shp' using driver `ESRI Shapefile'

NA Creating or updating layer places_merged failed.

NA Error: Write error.

buildings_merged = buildings_hospital$osm_points |>
  bind_rows(buildings_residential$osm_points, buildings_workship$osm_points) |>
  group_by(across(-geometry)) |> sf::st_cast("POINT") |>
  summarise(geometry = st_union(geometry), .groups = "drop")
sf::st_write(buildings_merged, "./03_Spatial_Data/POP/buildings_merged.shp", delete_dsn=T) 

NA Deleting source `./03_Spatial_Data/POP/buildings_merged.shp' failed
NA Writing layer `buildings_merged' to data source 
NA   `./03_Spatial_Data/POP/buildings_merged.shp' using driver `ESRI Shapefile'

NA Creating or updating layer buildings_merged failed.

NA Error: Write error.

# Visualize
tmap::tm_shape(aoi_country) + tmap::tm_borders(col="purple", lwd=1) +
  tmap::tm_shape(aoi_states) + tmap::tm_borders(col="purple", lwd=.5) +
#  tmap::tm_shape(buildings_merged) + tmap::tm_symbols(size=0.15,lwd=0.5,fill="white",col="purple") +
#  tmap::tm_add_legend(type="symbols", col="white", fill="purple", size=0.8, labels="Buildings") +
  tmap::tm_shape(places_merged) + tmap::tm_symbols(size=0.2,lwd=0.5,fill="white",col="blue") +
  tmap::tm_add_legend(type="symbols", col="white", fill="blue", size=0.8, labels="Towns & Villages") +
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
  tmap::tm_compass(color.dark="gray60",text.color="gray60",position=c("RIGHT", "top")) +
  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_layout(legend.position=c("left", "top"), legend.bg.color = "white") +
  tmap::tm_title("Risk Covariates Map: Built Environment", size=.8) + 
  tmap::tm_basemap("OpenStreetMap") 

  #tmap::tm_basemap("OpenTopoMap") 
  #tmap::tm_basemap("Esri.WorldImagery") 

Annualized deforestation

# Assign zones
zones_sf = counties |> sf::st_transform("EPSG:32629")
zones_sf$zone_id <- 1:nrow(zones_sf)
zones_sv <- terra::vect(zones_sf)

# Calculate zonal annualization by jurisdiction
forest_2024 = terra::rast("./data/BINARY/forest_2024.tif")
forest_loss_2014_2019 = terra::rast("./data/BINARY/forest_loss_2014_2019.tif")
forest_loss_2019_2024 = terra::rast("./data/BINARY/forest_loss_2019_2024.tif")
forest_loss_2014_2024 = terra::rast("./data/BINARY/forest_loss_2014_2024.tif")
pixel_area_ha <- 0.088914  # 29.80124 x 29.80124 m² converted to hectares

zonal_2014_2019 <- terra::extract(forest_loss_2014_2019,zones_sv,fun=sum,na.rm=T)
zonal_2019_2024 <- terra::extract(forest_loss_2019_2024,zones_sv,fun=sum,na.rm=T)
names(zonal_2014_2019) <- c("zone_id", "loss_2014_2019")
names(zonal_2019_2024) <- c("zone_id", "loss_2019_2024")
zonal_2014_2019$loss_2014_2019 <- zonal_2014_2019$loss_2014_2019 * pixel_area_ha
zonal_2019_2024$loss_2019_2024 <- zonal_2019_2024$loss_2019_2024 * pixel_area_ha

# Merge baseline years
zonal_stats <- merge(
  zonal_2014_2019,      # e.g. (zone_id, loss_2014_2019)
  zonal_2019_2024,      # e.g. (zone_id, loss_2019_2024)
  by = "zone_id",       # Common ID column
  all = TRUE            # Keep all rows if zones differ
)
# Annualize 10-year total & rejoin to sf object
zonal_stats$loss_10yr <- zonal_stats$loss_2014_2019 + zonal_stats$loss_2019_2024
zonal_stats$annual_loss_10yr <- zonal_stats$loss_10yr / 10
zones_sf <- merge(zones_sf, zonal_stats, by="zone_id", all.x=TRUE)
head(zones_sf[, c("zone_id", "loss_2014_2019", "loss_2019_2024", 
                  "loss_10yr", "annual_loss_10yr")])

# Derive 10-yr annualized raster
zones_sv <- terra::vect(zones_sf)
annual_loss_10yr_raster <- rasterize(
  zones_sv,                  # polygon SpatVector
  forest_loss_2014_2019,     # template raster for resolution/extent
  field = "annual_loss_10yr") # the column to rasterize
names(annual_loss_10yr_raster) <- ("annual_loss_10yr")
raster::writeRaster(annual_loss_10yr_raster,"./data/BINARY/annual_loss_2014_2024_zonal.tif",overwrite=T)

Distance covariates

# Derive distance-to-edge raster 
forest_2024 = terra::rast("./data/BINARY/forest_2024.tif")
forest_mask <- terra::classify(forest_2024, cbind(c(1,0), c(1,NA)))
nonForest_mask <- terra::classify(forest_2024, cbind(c(1,0), c(NA,1)))
writeRaster(forest_mask, "./data/BINARY/forest_mask.tif", overwrite=T)
writeRaster(nonForest_mask, "./data/BINARY/nonForest_mask.tif", overwrite=T)

forest_for_distance <- forest_mask
forest_for_distance_raster <- raster::raster(forest_for_distance) 
distance_to_edge <- distance(forest_for_distance_raster, filename="./data/BINARY/distance_to_edge_unix.tif")

# Derive distance-to-feature rasters
template <- forest_2024
roads_rast <- rasterize(vect(transport), template, field=1, background=0)
roads_rast <- subst(roads_rast, from=0, to=NA)
distance_to_roads <- distance(roads_rast)
writeRaster(distance_to_roads, "./data/LULC/distance_to_roads.tif", overwrite=T)

places_vect   <- sf::st_read("./data/BINARY/places.shp") |> terra::vect()
places_rast   <- rasterize(places_vect,template,field=1,background=0,touches=T)
places_rast   <- subst(places_rast, from=0, to=NA)
distance_to_places<- distance(places_rast)
writeRaster(distance_to_places, "./data/LULC/distance_to_places.tif", overwrite=T)

waterways_rast<-rasterize(vect(waterways_rast), template, field=1, background=0)
waterways_rast<-subst(waterways_rast, from=0, to=NA)
distance_to_waterways<-distance(waterways_rast)
writeRaster(distance_to_waterways, "./data/LULC/distance_to_waterways.tif", overwrite=T)

urban_rast    <- raster::raster(urban)
urban_rast    <- subst(urban_rast, from=0, to=NA)
distance_to_urban <- distance(urban_rast)
writeRaster(distance_to_urban, "./data/LULC/distance_to_urban.tif", overwrite=T)

# Normalize covariates for quicker computing
normalize <- function(x){
  (x-global(x,"min",na.rm=T))/
    (global(x,"max",na.rm=T) -
       global(x, "min",na.rm=T))}

population= raster::raster(population)
population= normalize_function(population)
slope     = raster::raster(slope)
slope     = normalize_function(slope)

distance_to_forest    = normalize_function(distance_to_forest)
distance_to_roads     = normalize_function(distance_to_roads)
distance_to_places    = normalize_function(distance_to_places)
distance_to_water     = normalize_function(distance_to_waterways)
distance_to_urban     = normalize_function(distance_to_urban)

# invert all risk-producing covariates
distance_to_edge_inv <- 1 - distance_to_forest
distance_to_roads_inv  <- 1 - distance_to_roads
distance_to_places_inv <- 1 - distance_to_places
distance_to_water_inv <- 1 -  distance_to_waterways
distance_to_urban_inv   <- 1 - distance_to_urban

writeRaster(distance_to_edge_inv, "./data/BINARY/distance_to_edge_inverted.tif", overwrite=T)
writeRaster(distance_to_roads_inv, "./data/BINARY/distance_to_roads_inverted.tif", overwrite=T)
writeRaster(distance_to_places_inv, "./data/BINARY/distance_to_places_inverted.tif", overwrite=T)
writeRaster(distance_to_water_inv, "./data/BINARY/distance_to_water_inverted.tif", overwrite=T)
writeRaster(distance_to_urban_inv, "./data/BINARY/distance_to_urban_inverted.tif", overwrite=T)
writeRaster(slope_norm, "./data/BINARY/slope_norm.tif", overwrite=T)

Deforestation risk

Two methods were explored for weighting variables and creating a generalized deforestation risk index. We could consider developing a spatial risk model using the spatstat package or logistic regression, as has been cited in recent Verra guides. In addition, some of the heavy lifting with input formatting and data wrangling has already been completed.

However, spatial modelling has tended to produce challenges when fitting such large country-wide covariates. Moreover, these kinds of spatialy driven models tend to require longer training procedures due to their intercept-based spatial kernels and slower resampling patterns.

Alternatively, we have drafted a tentative risk indexing approach based on a weighted sum of subjectively scored covariate effects. While each variable would still need a carefully reasoned score, this option offers a more streamlined method that is easier to adjust. We applied this risk index to inform a risk weighted allocation of the 10-year deforestation rate, first by multiplying the fraction of pixel risk by zonal forest loss, and second by factoring out annual zonal loss by multiplying by pixel risk values, as shown below

\[ \mathrm{AllocatedLoss}_{\mathrm{pixel}} = \left( \frac{\mathrm{risk}_{\mathrm{pixel}}}{\sum \mathrm{risk}_{\mathrm{zone}}} \right) \times \mathrm{annual\_loss\_10yr}_{\mathrm{zone}} \]

\[ \mathrm{allocated\_loss}_{\mathrm{pixel}} = \mathrm{risk}_{\mathrm{pixel}} \times \left( \frac{\mathrm{annual\_loss\_10yr}_{\mathrm{zone}}}{\sum \mathrm{risk}_{\mathrm{zone}}} \right) \]

Both formulas describe the same operation in different orders of multiplication: each pixel in a given zone Z receives a share of annual_loss_10yr_Z based on its proportional risk (the pixel’s risk relative to the sum of all pixel risks in that zone). This ensures that higher-risk pixels are allocated more deforestation, in line with the Verra guidance for an allocated deforestation risk map.

We intend to present both of these approaches for broader review and discussion in our upcoming meeting.

risk_index <- (0.2 * distance_to_edge) +
              (0.2 * distance_to_roads) +
              (0.2 * distance_to_places) +
              (0.1 * distance_to_urban) +
              (0.1 * distance_to_water) +
              (0.1 * slope)

# Re-normalize the index to between 0 and 1
rmin <- global(risk_index, "min", na.rm=TRUE)[1]
rmax <- global(risk_index, "max", na.rm=TRUE)[1]
risk_index_norm <- (risk_index - rmin) / (rmax - rmin)
writeRaster(risk_index_norm, "./data/BINARY/deforestation_risk_index.tif",overwrite=T)

# Returns a data.frame with columns: ID, risk_index_norm_sum
risk_sum <- extract(risk_index_norm, zones_sv, fun = sum, na.rm=TRUE)
colnames(risk_sum) <- c("zone_id","sum_risk")
zones_stats <- merge(zones_sf, risk_sum, by="zone_id", all.x=TRUE)
zones_rast <- rasterize(vect(zones_stats), risk_index_norm, field="zone_id")
zones_stats$loss_factor <- zones_stats$annual_loss_10yr / zones_stats$sum_risk
loss_factor_rast <- rasterize(vect(zones_stats), risk_index_norm, field="loss_factor")
allocated_loss <- risk_index_norm * loss_factor_rast
writeRaster(allocated_loss, "./data/BINARY/allocated_deforestation.tif", overwrite=T)

# Visualize
terra::plot(risk_index_norm, main="Deforestation risk map")
terra::plot(allocated_loss, main="Allocated deforestation map")

Appendix I: Reference Period Classifiers

To run these, you may change eval=F to eval=T at the top of chunk in the .Rmd or .R file saved in the OneDrive folder.

########################### 2019
# extract signatures
signatures_2019 = raster::extract(STACK_2019, samples ,df=T) # watch for data formats
samples_signatures_2019 <- dplyr::inner_join(signatures_2019, samples, by=c("ID"="id"))
samples_signatures_2019$geometry <- NULL # set geometry to NULL for model training

# training-test split, p=0.7 -> 70% split
partitioned_data_2019 <- stratified_partition(samples_signatures_2019,group_col="label", train_ratio=0.7)
trainData_2019 <- partitioned_data_2019$train
testData_2019 <- partitioned_data_2019$test
table(trainData_2019$label)
table(testData_2019$label)

# synthetic minority oversampling technique
trainData_2019<-performanceEstimation::smote(label ~ .,data=trainData_2019,perc.over=10, perc.under=100)
testData_2019<-performanceEstimation::smote(label ~ .,data=testData_2019,perc.over=10,perc.under=100)
# interpolate NAs with class-median-normalization (NAs -> missing cloud pixels)
trainData_2019 <- trainData_2019 |> group_by(label) |> mutate(across(where(is.numeric),
    ~ ifelse(is.na(.), median(., na.rm = TRUE), .))) |> ungroup()
testData_2019 <- testData_2019 |> group_by(label) |> mutate(across(where(is.numeric),
    ~ ifelse(is.na(.), median(., na.rm = TRUE), .))) |> ungroup()

# assign model variables
response  <- c("label")
predictors_2019 <- c(
  "NDVI_2019", "BLUE_2019", "GREEN_2019", "RED_2019", 
  "NIR08_2019", "SWIR16_2019", "SWIR22_2019", "DEM"
  )

# train classifier
rf_model_2019 <- caret::train(
  label~.,
  data = trainData_2019[, c(predictors_2019, "label")], # drop ID var
  trControl = cv_regime,
  method    = "rf", 
  metric    = 'Kappa', 
  ntree     = 500,
  tuneLength= 6,
  importance= T
  )

rf_test_2019 <- predict(rf_model_2019, testData_2019)
print(rf_model_2019) # cv results
confusionMatrix(rf_test_2019,testData_2019$label) # blind test results

index_feature_2019 <- createMultiFolds(trainData_2019$label, times=5) 
predictor_seq_2019 <-seq(from=1, to=length(predictors_2019),by=2)

subset_regime_2019 <- rfeControl(
  method="LGOCV",
  number = 10,
  verbose=FALSE,
  functions=rfFuncs,
  index=index_feature_2019
  )

rf_model_subset_2019 <- caret::rfe(
  label~.,
  data = trainData_2019[, c(predictors_2019, "label")], 
  sizes = predictor_seq_2019,
  metric = "Kappa",
  ntree=500,
  method="rf",
  rfeControl = subset_regime_2019
  )

rf_subset_test_2019 <- predict(rf_model_subset_2019,testData_2019)
print(rf_model_subset_2019)
confusionMatrix(rf_subset_test_2019$pred,testData_2019$label)

######################### 2024
# extract signatures
signatures_2024 = raster::extract(STACK_2024, samples ,df=T) # watch for data formats
samples_signatures_2024 <- dplyr::inner_join(signatures_2024, samples, by=c("ID"="id"))
samples_signatures_2024$geometry <- NULL # set geometry to NULL for model training

# training-test split, p=0.7 -> 70% split
partitioned_data_2024 <- stratified_partition(samples_signatures_2024,group_col="label", train_ratio=0.7)
trainData_2024 <- partitioned_data_2024$train
testData_2024 <- partitioned_data_2024$test
table(trainData_2024$label)
table(testData_2024$label)

# synthetic minority oversampling technique
trainData_2024<-performanceEstimation::smote(label ~ .,data=trainData_2024,perc.over=10, perc.under=100)
testData_2024<-performanceEstimation::smote(label ~ .,data=testData_2024,perc.over=10,perc.under=100)
# interpolate missing cloud pixels with class-median-normalization
trainData_2024 <- trainData_2024 |> group_by(label) |> mutate(across(where(is.numeric),
    ~ ifelse(is.na(.), median(., na.rm = TRUE), .))) |> ungroup()
testData_2024 <- testData_2024 |> group_by(label) |> mutate(across(where(is.numeric),
    ~ ifelse(is.na(.), median(., na.rm = TRUE), .))) |> ungroup()

# interpolate NAs with class-median-normalization (NAs -> missing cloud pixels)
trainData_2024 <- trainData_2024 |> group_by(label) |> mutate(across(where(is.numeric),
    ~ ifelse(is.na(.), median(., na.rm = TRUE), .))) |> ungroup()
testData_2024 <- testData_2024 |> group_by(label) |> mutate(across(where(is.numeric),
    ~ ifelse(is.na(.), median(., na.rm = TRUE), .))) |> ungroup()
water_2014 <- trainData_2014[trainData_2014$label == "Water", ]
water_ids <- water_2014$ID
water_2024 <- samples_signatures_2024[samples_signatures_2024$ID %in% water_ids, ]
trainData_2024 <- rbind(trainData_2024, water_2024)
table(trainData_2024$label)

# assign model variables
response  <- c("label")
predictors_2024 <- c(
  "NDVI_2024", "BLUE_2024", "GREEN_2024", "RED_2024", 
  "NIR08_2024", "SWIR16_2024", "SWIR22_2024", "DEM"
  )

# train classifier
rf_model_2024 <- caret::train(
  label~.,
  data = trainData_2024[, c(predictors_2024, "label")], # drop ID var
  trControl = cv_regime,
  method    = "rf", 
  metric    = 'Kappa', 
  ntree     = 500,
  tuneLength= 6,
  importance= T
  )

rf_test_2024 <- predict(rf_model_2024, testData_2024)
print(rf_model_2024) # cv results
confusionMatrix(rf_test_2024,testData_2024$label) # blind test results

index_feature_2024 <- createMultiFolds(trainData_2024$label, times=5) 
predictor_seq_2024 <-seq(from=1, to=length(predictors_2024),by=2)

subset_regime_2024 <- rfeControl(
  method="LGOCV",
  number = 10,
  verbose=FALSE,
  functions=rfFuncs,
  index=index_feature_2024
  )

rf_model_subset_2024 <- caret::rfe(
  label~.,
  data = trainData_2024[, c(predictors_2024, "label")], 
  sizes = predictor_seq_2024,
  metric = "Kappa",
  ntree=500,
  method="rf",
  rfeControl = subset_regime_2024
  )

rf_subset_test_2024 <- predict(rf_model_subset_2024,testData_2024)
print(rf_model_subset_2024)
confusionMatrix(rf_subset_test_2024$pred,testData_2024$label)

LULC_LIBERIA_2014 <- raster::predict(STACK_2014,rf_model_2014, na.rm=TRUE) |> raster::raster()
raster::writeRaster(LULC_LIBERIA_2014,"./data/LULC/LULC_LIBERIA_2014-01-04.tif",
  format = "GTiff", overwrite = T)

LULC_LIBERIA_2019 <- raster::predict(STACK_2019,rf_model_2019, na.rm=TRUE) |> raster::raster() 
raster::writeRaster(LULC_LIBERIA_2019,"./data/LULC/LULC_LIBERIA_2019-01-02.tif",
  format = "GTiff",overwrite = T)

LULC_LIBERIA_2024 <- raster::predict(STACK_2024,rf_model_2014, na.rm=TRUE) |> raster::raster() 
raster::writeRaster(LULC_LIBERIA_2024,"./data/LULC/LULC_LIBERIA_2024-01-16x.tif",
  format = "GTiff",overwrite = T)

Appendix II: Runtime snapshot

devtools::session_info()

LULC_PROJECT_2015=terra::rast("./data/LULC/LULC_PROJECT_2015.tif")
LULC_PROJECT_2019=terra::rast("./data/LULC/LULC_PROJECT_2019.tif")
LULC_PROJECT_2024=terra::rast("./data/LULC/LULC_PROJECT_2024.tif")
LULC_PROJECT_2015=terra::project(LULC_PROJECT_2015, "EPSG:32629")
LULC_PROJECT_2024=terra::project(LULC_PROJECT_2024, "EPSG:32629")
LULC_PROJECT_2019=terra::project(LULC_PROJECT_2019, "EPSG:32629")

voi = sf::st_transform(aoi, 32629) |> terra::vect() 
LULC_PROJECT_2015 = terra::crop(LULC_PROJECT_2015, voi, mask=T)
LULC_PROJECT_2019 = terra::crop(LULC_PROJECT_2019, voi, mask=T)
LULC_PROJECT_2024 = terra::crop(LULC_PROJECT_2024, voi, mask=T)

code_dict <- data.frame(
  id = c(0, 1, 2, 3, 4, 5, 6, 7),
  label = c("Water", "Forest", "Grassland", 
            "Wetland", "Croplan", "Shrubland", 
            "Urban", "Bareground"))

levels(LULC_PROJECT_2015) <- code_dict
levels(LULC_PROJECT_2019) <- code_dict
levels(LULC_PROJECT_2024) <- code_dict

terra::plot(LULC_PROJECT_2015, main = "LULC 2015")
terra::plot(LULC_PROJECT_2015, main = "LULC 2019")
terra::plot(LULC_PROJECT_2015, main = "LULC 2024")

raster_res <- terra::res(LULC_PROJECT_2015) # 9.933065 x 9.933065 m2
pixel_area_ha <- (raster_res[1] * raster_res[2]) / 10000  # Convert m² to hectares
freq_2015 <- as.data.frame(terra::freq(LULC_PROJECT_2015))
freq_2015$area_ha <- freq_2015$count * pixel_area_ha
freq_2015$percentage <- (freq_2015$area_ha / sum(freq_2015$area_ha)) * 100
freq_2015$year <- 2015

freq_2019 <- as.data.frame(terra::freq(LULC_PROJECT_2019))
freq_2019$area_ha <- freq_2019$count * pixel_area_ha
freq_2019$percentage <- (freq_2019$area_ha / sum(freq_2019$area_ha)) * 100
freq_2019$year <- 2019

freq_2024 <- as.data.frame(terra::freq(LULC_PROJECT_2024))
freq_2024$area_ha <- freq_2024$count * pixel_area_ha
freq_2024$percentage <- (freq_2024$area_ha / sum(freq_2024$area_ha)) * 100
freq_2024$year <- 2024

land_cover_summary=bind_rows(freq_2015, freq_2019, freq_2024)
land_cover_summary=merge(land_cover_summary,code_dict, by.x="value",by.y="id",all.x=T)

land_cover_summary_wide <- land_cover_summary %>%
  select(label, year, area_ha, percentage) %>%
  pivot_wider(names_from = year, values_from = c(area_ha, percentage))
print(land_cover_summary_wide)

unique_values_2015 <- terra::unique(LULC_PROJECT_2015)
unique_values_2019 <- terra::unique(LULC_PROJECT_2019)
unique_values_2024 <- terra::unique(LULC_PROJECT_2024)

print(unique_values_2015)
print(unique_values_2019)
print(unique_values_2024)

print(terra::is.factor(LULC_PROJECT_2015))
print(terra::levels(LULC_PROJECT_2015))
print(terra::levels(LULC_PROJECT_2019))
print(terra::levels(LULC_PROJECT_2024))


levels(LULC_PROJECT_2015) <- data.frame(ID = code_dict$id, label = code_dict$label)
levels(LULC_PROJECT_2019) <- data.frame(ID = code_dict$id, label = code_dict$label)
levels(LULC_PROJECT_2024) <- data.frame(ID = code_dict$id, label = code_dict$label)
freq_2015 <- as.data.frame(terra::freq(LULC_PROJECT_2015))
freq_2019 <- as.data.frame(terra::freq(LULC_PROJECT_2019))
freq_2024 <- as.data.frame(terra::freq(LULC_PROJECT_2024))

print(freq_2015)
print(freq_2019)
print(freq_2024)

print(LULC_PROJECT_2015)
print(LULC_PROJECT_2019)
print(LULC_PROJECT_2024)

plot(LULC_PROJECT_2015, main = "LULC 2015")
plot(LULC_PROJECT_2019, main = "LULC 2019")
plot(LULC_PROJECT_2024, main = "LULC 2024")

References

Government of Liberia. 2019. “Liberia’s Forest Reference Emission Level Submission to the UNFCCC.”